library("rstudioapi")     
##Set working directory
setwd(dirname(getActiveDocumentContext()$path))
getwd()
library("readxl") 
require(miceadds)
library(data.table)
library(tidyverse)
library(DescTools)
library(htmlTable)
library(stargazer)
library(xlsx)
library(estimatr)

# Load final Data file: 
load("FinalFinal20062023Full1-6.RData")
load("FinalFinalSP20062023Full1-6.RData")

library(dplyr)
final_finalV2 = mutate(final_finalV2, 
                       Agec = (Age - 30)/10, # per 10 year age
                       Incomec = (Income - 100)/50, # per 50 income
                       clinic_distancec = (clinic_distance - 5)/5, # per 5
                       SNMetricc = SNMetric/10)
final_finalSP = mutate(final_finalSP, 
                       Agec = (Age - 30)/10, # per 10 year age
                       Incomec = (Income - 100)/50, # per 50 income
                       clinic_distancec = (clinic_distance - 5)/5, # per 5
                       SNMetricc = SNMetric/10)

final_finalV2$Village_Population_Size = relevel(as.factor(final_finalV2$Village_Population_Size), ref = "Small")

L_8dV2 <- glm.cluster(vaccine_reported_combo ~ Ncash_catV2*Village_Population_Size + Gender + Agec + 
                      Education + Incomec + Employed + SNMetricc + clinic_distancec+
                      dist2 + dist3 + dist4 + dist5 + dist6, 
                    data = final_finalV2, family=binomial,
                    cluster = final_finalV2$Q123) 
exp(cbind("Odds ratio" = coef(L_8dV2), confint.default(L_8dV2, level = 0.95)))


L_9dV2 <- glm.cluster(ActVacApril  ~ Ncash_catV2*Village_Population_Size + Gender + Agec + 
                      Education + Incomec + Employed + SNMetricc + clinic_distancec+
                      dist2 + dist3 + dist4 + dist5 + dist6, 
                    data = final_finalV2, family=binomial,
                    cluster = final_finalV2$Q123) #### individual treatment
exp(cbind("Odds ratio" = coef(L_9dV2), confint.default(L_8dV2, level = 0.95)))


tab_8d1 = exp(cbind("Odds ratio" = coef(L_8dV2), confint.default(L_8dV2, level = 0.95)))
tab_9d1 = exp(cbind("Odds ratio" = coef(L_9dV2), confint.default(L_9dV2, level = 0.95)))
rownames(tab_8d1)[c(1:23)] =  c("Intercept", "Health", "Health_Placebo", "HighCash", "HighCash_Placebo", "LowCash", "LowCash_Placebo","Village_Population_SizeLarge","Village_Population_SizeMid", "Male", "Age", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food", "Employed", "Social Media", "Access","Discrict 2", "Discrict 3", "Discrict 4", "Discrict 5","Discrict 6")
rownames(tab_9d1)[c(1:23)] = c("Intercept", "Health", "Health_Placebo", "HighCash", "HighCash_Placebo", "LowCash", "LowCash_Placebo", "Village_Population_SizeLarge","Village_Population_SizeMid" ,"Male", "Age", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food", "Employed", "Social Media", "Access","Discrict 2", "Discrict 3", "Discrict 4", "Discrict 5","Discrict 6")

rownames(tab_8d1)[c(24:35)] = c("CDC_SizeLarge","CDC_placebo_SizeLarge",
                                "HighCash_SizeLarge", "HighCash_placebo_SizeLarge",
                                "LowCash_SizeLarge", "LowCash_Placebo_SizeLarge",
                                "CDC_SizeMid", "CDC_Placebo_SizeMid", 
                                "HighCash_SizeMid", "HighCash_placebo_SizeMid",
                                "LowCash_SizeMid", "LowCash_Placebo_SizeMid")
rownames(tab_9d1)[c(24:35)] = c("CDC_SizeLarge","CDC_placebo_SizeLarge",
                                "HighCash_SizeLarge", "HighCash_placebo_SizeLarge",
                                "LowCash_SizeLarge", "LowCash_Placebo_SizeLarge",
                                "CDC_SizeMid", "CDC_Placebo_SizeMid", 
                                "HighCash_SizeMid", "HighCash_placebo_SizeMid",
                                "LowCash_SizeMid", "LowCash_Placebo_SizeMid")

convert_fun = function(data, name_mod){
  new_vec = c()
  for (i in c(1:dim(data)[1])){
    new_vec[i] = as.character(paste(format(round(data[i,1],2), nsmall = 2),paste("(",paste(format(round(data[i,2],2), nsmall = 2), format(round(data[i,3],2), nsmall = 2), sep = ", "),")", sep ="")))
  }
  new_vec = data.frame(new_vec)
  rownames(new_vec) = rownames(data)
  colnames(new_vec) = name_mod
  return(new_vec)
}


final_finalSP$Village_Population_Size = relevel(as.factor(final_finalSP$Village_Population_Size), ref = "Small")


L_7csV2 <- glm.cluster(vaccine_reported ~ Vcash_catV2*Village_Population_Size + Gender + Agec + 
                       Education + Incomec + Employed + SNMetricc + clinic_distancec+
                       dist2 + dist3 + dist4 + dist5 + dist6, 
                     data = final_finalSP, family=binomial, 
                     cluster = final_finalSP$Q123) 

L_8csV2 <- glm.cluster(ActVacApril ~ Vcash_catV2*Village_Population_Size + Gender + Agec + 
                       Education + Incomec + Employed + SNMetricc + clinic_distancec+
                       dist2 + dist3 + dist4 + dist5 + dist6, 
                     data = final_finalSP, family=binomial, 
                     cluster = final_finalSP$Q123) 


tab_7cs = exp(cbind("Odds ratio" = coef(L_7csV2), confint.default(L_7csV2, level = 0.95)))
tab_8cs = exp(cbind("Odds ratio" = coef(L_8csV2), confint.default(L_8csV2, level = 0.95)))

rownames(tab_7cs)[c(1:20)] = c("Intercept", "VCash_Health","VCash_HighCash", "VCash_LowCash", "Village_Population_SizeLarge","Village_Population_SizeMid" , "Male", "Age", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food", "Employed", "Social Media", "Access", "Discrict 2", "Discrict 3", "Discrict 4", "Discrict 5","Discrict 6")
rownames(tab_8cs)[c(1:20)] = c("Intercept", "VCash_Health","VCash_HighCash", "VCash_LowCash", "Village_Population_SizeLarge","Village_Population_SizeMid" , "Male", "Age", "High-Educate", "Low-Educate", "Medium-Educate", "Mean-Food", "Employed", "Social Media", "Access", "Discrict 2", "Discrict 3", "Discrict 4", "Discrict 5","Discrict 6")
rownames(tab_7cs)[c(21:26)] = c("CDC_SizeLarge", "HighCash_SizeLarge", "LowCash_SizeLarge",
                                "CDC_SizeMid", "HighCash_SizeMid", "LowCash_SizeMid")
rownames(tab_8cs)[c(21:26)] = c("CDC_SizeLarge", "HighCash_SizeLarge", "LowCash_SizeLarge",
                                "CDC_SizeMid", "HighCash_SizeMid", "LowCash_SizeMid")

col_8d1 = convert_fun(tab_8d1, "Model1")
col_9d1 = convert_fun(tab_9d1, "Model2")
col_7cs = convert_fun(tab_7cs, "ModelSP1")
col_8cs = convert_fun(tab_8cs, "ModelSP2")

tsp1 = merge(col_8d1, col_9d1, by = 0, all = TRUE);rownames(tsp1) = tsp1[,1];tsp1[,1] = NULL
tsp2 = merge(col_7cs, col_8cs, by = 0, all = TRUE);rownames(tsp2) = tsp2[,1];tsp2[,1] = NULL
tbsp1 = merge(tsp1, tsp2, by = 0, all = TRUE) ;rownames(tbsp1) = tbsp1[,1];tbsp1[,1] = NULL

tbsp2 = tbsp1[c(3,4,5,6,13,14,16,17,18,19,20,21,24,25,26,27,28,29,34,35,36,37,38,22),]
tbsp2[25,] = rep("Yes", 4)
rownames(tbsp2)[25] = "Covariates"
tbsp2[26,] = rep("Yes",4)
rownames(tbsp2)[26] = "Fixed Effects"


no.obs = c((L_8dV2$glm_res$df.null), (L_9dV2$glm_res$df.null), (L_7csV2$glm_res$df.null), (L_8csV2$glm_res$df.null))
aic_mod = round(c(L_8dV2$glm_res$aic, L_9dV2$glm_res$aic, L_7csV2$glm_res$aic, L_8csV2$glm_res$aic),2)
log_lik = round(c(as.vector(logLik(L_8dV2$glm_res)), as.vector(logLik(L_9dV2$glm_res)), as.vector(logLik(L_7csV2$glm_res)), as.vector(logLik(L_8csV2$glm_res))),2)
mod_spec = data.frame(rbind(no.obs, aic_mod, log_lik))
rownames(mod_spec) = c("Observations", "Akaike Inf. Crit.", "Log Likelihood")
colnames(mod_spec) = colnames(tbsp2)

tbsp3 = rbind(tbsp2, mod_spec)

require(xtable)
xtable(tbsp3)

# Table S5: 
stargazer(tbsp3, type = "text", summary = FALSE, out = "SM_TableS5.txt")
stargazer(tbsp3, type = "latex", summary = FALSE, out = "SM_TableS5.txt")

